# Replication code for Ryan Brutger and Amy Pond's PSRM manuscript (Last updated March 10, 2025)
# This r script generates the results for the manuscript and appendix for the survey results

# Code was run on R version R version 3.6.1 (2019-07-05) -- "Action of the Toes"
# Using MacBook Pro using macOS Ventura, with 2.8 GHz Quad-Core Intel Core i7

# To generate the marginal effects plats, you must first run the companion R script
# "interaction_plots.R" to generate functions created by Anton Strezhnev (06/17/2013)

################# 1. Load libraries, load and clean data ##############

# First, set your working directory using the setwd() command:
setwd()

#install.packages("pacman")
pacman::p_load(dplyr, ggplot2, quanteda, plotly, foreign, interplot, plm, reshape2, lda,
               countrycode, sandwich, lmtest, MASS, RColorBrewer, states, mice, VIM, 
               margins, clusterSEs, optimx, coefplot, systemfit, randomizr, estimatr, ri2,
               stargazer, fastDummies, cowplot, egg, cobalt, readr, xtable, sjPlot, sjmisc) # load packages


#### Data Pre-Processing #### 
data <- read_csv("Antitrust_PSRM_Data_Module_1.csv", col_types = cols (.default = "i", Treat1 = "c", DomFor = "c", DomFor2 = "c"))


# Support in control and full sample mentioned in the text
mean(data$AT_StrengthPos[data$Control==1], na.rm=TRUE) #59%

# Among Democrats
mean(data$AT_StrengthPos[data$Control==1 & data$dem==1], na.rm=TRUE) #67%

# Among Republicans
mean(data$AT_StrengthPos[data$Control==1 & data$rep==1], na.rm=TRUE) #59%


# Differences between Rep and Dems
DemRepData <- subset(data[data$party<3, ]) # subset to only those who identify as Democrat or Republican
DemRepStrengthen <- lm(AT_StrengthPos ~ dem, data = DemRepData[DemRepData$Control==1, ]) # Test partisan difference in the control condition
summary(DemRepStrengthen) # 8.27, p=0.032 reported in the text

# Models for Table 1 of manuscript
AT.Strength <- lm(AT_Strength ~ CompPrices + CompSmall + Punish, data = data) #Main Effects 
AT.Strength.Cont <- lm(AT_Strength ~ CompPrices + CompSmall +  Punish + rep + college + male + hh_income, data = data) #Main Effects with controls
AT.Strength.int <- lm(AT_Strength ~ CompPrices*rep + CompSmall*rep +  Punish*rep, data = data[data$party<3, ])  #Main Effects  with Republican interaction, Democrats are baseline
AT.Strength.int.c <- lm(AT_Strength ~ CompPrices*rep + CompSmall*rep +  Punish*rep + college + male + hh_income, data = data[data$party<3, ])#Main Effects with Republican interaction and controls, Democrats are baseline

# Create Table for all module 1 results
# Note - signs for significance were manually updated so "cross" is for p<0.1, * for p<0.05, and ** for p<0.01
stargazer(AT.Strength, AT.Strength.Cont, AT.Strength.int, AT.Strength.int.c, omit.stat = c("rsq", "adj.rsq", "ser", "f"))

#
# For marginal effects, first run "interaction_plots.R" to generate functions created by Anton Strezhnev 06/17/2013
#

#Create Figure 2 of manuscript: marginal effects of Punish Companies treatment for Democrats and Republicans
# manually set ylim for plot to ylim(-0.5, 0.01)
# export with height = 3.5, width = 4
par(mar=c(2,5,3,1))
interaction_plot_binary(AT.Strength.int, effect="Punish", moderator="rep", 
                        interaction="rep:Punish", factor_labels=c("Democrat","Republican"), 
                        xlabel="", ylabel="Effect of Punish Companies on \n Support for Strengthening Antitrust Law", 
                        title="")


# Create Table A3 of the appendix: Results with Domestic Wording Only
# Note - signs for significance were manually updated so "cross" is for p<0.1, * for p<0.05, and ** for p<0.01
AT.Strength.int.for <- lm(AT_Strength ~ CompPrices + CompSmall  + Punish, data = data[data$Foreign==0, ])
AT.Strength.int.for.c <- lm(AT_Strength ~ CompPrices + CompSmall + Punish + rep + college + male + hh_income, data = data[data$Foreign==0, ])
stargazer(AT.Strength.int.for, AT.Strength.int.for.c, omit.stat = c("rsq", "adj.rsq", "ser", "f"), column.labels = c("Strengthen Laws", "Economy"))

# Create Table A4 of the appendix: Tech Interactions
# Note - signs for significance were manually updated so "cross" is for p<0.1, * for p<0.05, and ** for p<0.01
AT.Strength.int.tech <- lm(AT_Strength ~ CompPrices*Tech + CompSmall*Tech +  Punish*Tech, data = data)
AT.Strength.int.c.tech <- lm(AT_Strength ~ CompPrices*Tech + CompSmall*Tech +  Punish*Tech + college + male + hh_income, data = data)
stargazer(AT.Strength.int.tech, AT.Strength.int.c.tech, omit.stat = c("rsq", "adj.rsq", "ser", "f"), column.labels = c("Strengthen Laws", "Economy"))

# Create Table A5 of the appendix: College Education Interactions for appendix
# Note - signs for significance were manually updated so "cross" is for p<0.1, * for p<0.05, and ** for p<0.01
AT.Strength.int.edu <- lm(AT_Strength ~ CompPrices*college + CompSmall*college +  Punish*college, data = data)
AT.Strength.int.c.edu <- lm(AT_Strength ~ CompPrices*college + CompSmall*college +  Punish*college + male + hh_income, data = data)
stargazer(AT.Strength.int.edu, AT.Strength.int.c.edu, omit.stat = c("rsq", "adj.rsq", "ser", "f"), column.labels = c("Strengthen Laws", "Economy"))

# Create Table A6 of the appendix: Political Knowledge Interactions 
# Note - signs for significance were manually updated so "cross" is for p<0.1, * for p<0.05, and ** for p<0.01
# Measures of political knowledge 
data$PK.brit <- ifelse(data$pol_knwldge_brit==3, 1, 0)
data$PK.house <- ifelse(data$pol_knwldge_house==1, 1, 0)
data$PK.tot <- data$PK.brit + data$PK.house # create consolidated measure of political knowledge
data$PK.high <- ifelse(data$PK.tot>1, 1, 0) # create high/low indicator for political knowledge

AT.Strength.int.pk <- lm(AT_Strength ~ CompPrices*PK.high + CompSmall*PK.high +  Punish*PK.high, data = data)
AT.Strength.int.c.pk <- lm(AT_Strength ~ CompPrices*PK.high + CompSmall*PK.high +  Punish*PK.high + college + male + hh_income, data = data)
stargazer(AT.Strength.int.pk, AT.Strength.int.c.pk, omit.stat = c("rsq", "adj.rsq", "ser", "f"), column.labels = c("Strengthen Laws", "Economy"))

# For section A.1 of the appendix: Module 1 DEMOGRAPHICS 
data$college <- ifelse(data$edu>3, 1, 0) #has a 2 or 4 yr degree or higher

data$age_18to24 <- ifelse(data$age >=18 & data$age <25,1,0)
data$age_25to39 <- ifelse(data$age >=25 & data$age <40, 1,0)
data$age_40to59 <- ifelse(data$age >= 40 & data$age <60, 1,0)
data$age_60plus <- ifelse(data$age >=60, 1,0)

data$less50k <- ifelse(data$hh_income<=2,1,0)
data$btw50_100k <- ifelse(data$hh_income>2 & data$hh_income<=4, 1,0)
data$btw100_150k <- ifelse(data$hh_income>4 & data$hh_income <=6, 1,0)
data$above150k <- ifelse(data$hh_income>=7,1,0)

demo_names <- c("Age 18 to 24",
                "Age 25 to 39",
                "Age 40 to 59",
                "Age >50",
                "Female",
                "Household income $0 to $50,000",
                "Household income $50,001 to $100,000",
                "Household income $100,001 to $150,000",
                "Household income >$150,000",
                "Attended college")
demo_vars <- list(data$age_18to24,
                  data$age_25to39,
                  data$age_40to59,
                  data$age_60plus,
                  data$gender1, # 1 = male, 2 = female
                  data$less50k,
                  data$btw50_100k,
                  data$btw100_150k,
                  data$above150k,
                  data$college)


pop_pct <- sprintf("%.3f",
                   round(c(# age breakdown from US Census Bureau 2019 Table 1: https://www.census.gov/data/tables/2019/demo/age-and-sex/2019-age-sex-composition.html (I ignore those <18)
                     0.1321,
                     0.2660,
                     0.3251,
                     0.2929,
                     .5097, # sex breakdown from US Census Bureau, 2019
                     .3707, #https://www.census.gov/data/tables/time-series/demo/income-poverty/cps-hinc/hinc-01.html
                     .2884,
                     .1555,
                     .1854,
                     0.6108 
                   ),digits=3)) # https://www.census.gov/content/census/en/data/tables/2019/demo/educational-attainment/cps-detailed-tables.html

prop.function <- function(data){
  prop <- as.numeric(prop.table(table(data))[2])
  prop <- sprintf("%.3f",round(prop,digits=3))
  return(prop)
}             

props <- lapply(demo_vars, prop.function)

demo_props <- as.data.frame(cbind(demo_names, props, pop_pct))
colnames(demo_props) <- c("Demographic", "Portion of US Sample",  "U.S. Population")
print(demo_props, row.names=F)

print(xtable(demo_props,
             caption = "Study demographics. U.S. population information on age, sex, income, and education are from the Census Bureau and are for 2019. Partisan identification is from Pew and covers registered voters for 2018/19."),
      include.rownames=F)

##
## Module 2 analysis
##

data2 <- read_csv("Antitrust_PSRM_Data_Module_2.csv", col_types = cols (.default = "i", DomFor2 = "c"))

# Create dichotomous dependent variable do those who think antitrust laws should be "somehwat" or "dramatically" strengthened
data2$F_StrengthPos <- ifelse(data2$F_Strength>3, 1, 0)
data2$F_StrengthPos <- data2$F_StrengthPos*100

#Effects of  Equality and Democracy treatments 
F.Strength  <- lm(F_Strength ~ F_Equality_Pro + F_Equity_Pro + F_Dem + F_Equality_Con + F_Equity_Con, data = data2)
F.Strength.c  <- lm(F_Strength ~  F_Equality_Pro + F_Equity_Pro + F_Dem + F_Equality_Con + F_Equity_Con + rep + college + male + hh_income, data = data2)
# Republican Equality and Democracy Interactions
F.Strength.rep.int  <- lm(F_Strength ~ F_Equality_Pro*rep + F_Equity_Pro*rep + F_Dem*rep + F_Equality_Con*rep + F_Equity_Con*rep , data = data2[data2$party<3, ])
F.Strength.c.rep.int  <- lm(F_Strength ~ F_Equality_Pro*rep + F_Equity_Pro*rep + F_Dem*rep + F_Equality_Con*rep + F_Equity_Con*rep  + college + male + hh_income, data = data2[data2$party<3, ])

#Create Table 2 of the manuscript 
# Note - the Negative-Equality and Negative-Equity treatment effects are ommited from Table 2 of the manuscript, but are displayed in Tables A1 and A2 of the appendix
# Note - signs for significance were manually updated so "cross" is for p<0.1, * for p<0.05, and ** for p<0.01
stargazer(F.Strength, F.Strength.c, F.Strength.rep.int, F.Strength.c.rep.int,  omit.stat = c("rsq", "adj.rsq", "ser", "f"), column.labels = c(""))

#Create Table A1 of the appendix
# Note - signs for significance were manually updated so "cross" is for p<0.1, * for p<0.05, and ** for p<0.01
stargazer(F.Strength, F.Strength.c,  omit.stat = c("rsq", "adj.rsq", "ser", "f"), column.labels = c(""))
#Create Table A2 of the appendix
# Note - signs for significance were manually updated so "cross" is for p<0.1, * for p<0.05, and ** for p<0.01
stargazer(F.Strength.rep.int, F.Strength.c.rep.int,  omit.stat = c("rsq", "adj.rsq", "ser", "f"), column.labels = c(""))


#Effects in percentages for in-text discussion
F.Strength.Pos  <- lm(F_StrengthPos ~ F_Equality_Pro + F_Equity_Pro + F_Dem + F_Equality_Con + F_Equity_Con, data = data2)
summary(F.Strength.Pos)# equality=4.9%, equity=5.3%


#Create Figure 3 of manuscript: marginal effects of Democracy treatment for Democrats and Republicans
# manually set ylim for plot to ylim(-0.25, 0.3)
# export with height = 3.5, width = 4
par(mar=c(2,5,3,1))
interaction_plot_binary(F.Strength.rep.int, effect="F_Dem", moderator="rep", 
                        interaction="rep:F_Dem", factor_labels=c("Democrat","Republican"), 
                        xlabel="", ylabel="Effect of Democracy Treatment on \n Support for Strengthening Antitrust Law", 
                        title="")


# For section A.1 of the appendix: Module 2 DEMOGRAPHICS 
data2$college <- ifelse(data2$edu>3, 1, 0) #has a 2 or 4 yr degree or higher

data2$age_18to24 <- ifelse(data2$age >=18 & data2$age <25,1,0)
data2$age_25to39 <- ifelse(data2$age >=25 & data2$age <40, 1,0)
data2$age_40to59 <- ifelse(data2$age >= 40 & data2$age <60, 1,0)
data2$age_60plus <- ifelse(data2$age >=60, 1,0)

data2$less50k <- ifelse(data2$hh_income<=2,1,0)
data2$btw50_100k <- ifelse(data2$hh_income>2 & data2$hh_income<=4, 1,0)
data2$btw100_150k <- ifelse(data2$hh_income>4 & data2$hh_income <=6, 1,0)
data2$above150k <- ifelse(data2$hh_income>=7,1,0)

demo_names <- c("Age 18 to 24",
                "Age 25 to 39",
                "Age 40 to 59",
                "Age >50",
                "Female",
                "Household income $0 to $50,000",
                "Household income $50,001 to $100,000",
                "Household income $100,001 to $150,000",
                "Household income >$150,000",
                "Attended college")
demo_vars <- list(data2$age_18to24,
                  data2$age_25to39,
                  data2$age_40to59,
                  data2$age_60plus,
                  data2$gender1, # 1 = male, 2 = female
                  data2$less50k,
                  data2$btw50_100k,
                  data2$btw100_150k,
                  data2$above150k,
                  data2$college)


pop_pct <- sprintf("%.3f",
                   round(c(# age breakdown from US Census Bureau 2019 Table 1: https://www.census.gov/data2/tables/2019/demo/age-and-sex/2019-age-sex-composition.html (I ignore those <18)
                     0.1321,
                     0.2660,
                     0.3251,
                     0.2929,
                     .5097, # sex breakdown from US Census Bureau, 2019
                     .3707, #https://www.census.gov/data2/tables/time-series/demo/income-poverty/cps-hinc/hinc-01.html
                     .2884,
                     .1555,
                     .1854,
                     0.6108 
                   ),digits=3)) # https://www.census.gov/content/census/en/data2/tables/2019/demo/educational-attainment/cps-detailed-tables.html

prop.function <- function(data2){
  prop <- as.numeric(prop.table(table(data2))[2])
  prop <- sprintf("%.3f",round(prop,digits=3))
  return(prop)
}             

props <- lapply(demo_vars, prop.function)


demo_props <- as.data.frame(cbind(demo_names, props, pop_pct))
colnames(demo_props) <- c("Demographic", "Portion of US Sample",  "U.S. Population")
print(demo_props, row.names=F)

print(xtable(demo_props,
             caption = "Study demographics. U.S. population information on age, sex, income, and education are from the Census Bureau and are for 2019. Partisan identification is from Pew and covers registered voters for 2018/19."),
      include.rownames=F)

